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(57) ABSTRACT 


.An apparatus, computer program product and method of ana- 
lyzing non-stationary time varying phenomena. A represen- 
tation of a non-stationary time varying phenomenon is recur- 
sively sifted using Empirical Mode Decomposition (EMD) to 
extract intrinsic mode functions (IMFs). The representation is 
filtered to extract intrinsic trends by combining a number of 
IMFs. The intrinsic trend is inherent in the data and identifies 
an IMF indicating the variability of the phenomena. The trend 
also may be used to detrend the data. 
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ANALYZING NONSTATIONARY FINANCIAL 
TIME SERIES VIA HILBERT-HUANG 
TRANSFORM (HHT) 

CROSS REFERENCE TO RELATED 5 

APPLICATIONS 

The present application is a continuation of U.S. Provi- 
sional Patent Application Ser. No. 60/510,678 filed Oct. 9, 
2003 and related to U.S. Pat. No. 6,782,124 entitled “Three to 
Dimensional Empirical Mode Decomposition Analysis 
Apparatus Method and Article of Manufacture” to Per Glo- 
ersen, assigned to the assignee of the present invention. 

ORIGIN OF THE INVENTION 15 

The invention described herein was made by an employee 
of the United States Government, and may be manufactured 
and used by or for the Government for governmental pur- 
poses without the payment of any royalties thereon or there- 20 
for. 

BACKGROUND OF THE INVENTION 

1 . Field of the Invention 25 

The present invention is generally related to a data analysis 

method, apparatus and article of manufacture and more par- 
ticularly to an apparatus, article of manufacture and analysis 
method for determining and extracting trends from nonsta- 
tionary time varying data. 30 

Although the present invention finds utility in determining 
and extracting trends from non stationary time varying data 
sequences such as financial and climatological data, it is 
understood that the present invention has application to any 
non stationary time varying data representative of real world 35 
phenomenon. For example, the real world phenomenon to 
which the invention finds utility may include any of a wide 
variety of real world phenomena such as population growth, 
traffic flow and, non stationary time varying data representa- 
tive of processes including electrical, mechanical, biological, 40 
chemical, optical, geophysical or other process(es) that may 
be analyzed and thereby more fully understood by applying 
the invention thereto. 

2. Background Description 

Very often when large volumes of data are collected about 45 
real word processes, one of the first components analyzed is 
an apparently aperiodic component or, known as the data 
trend, e.g., “the markets are trending up.” For a rigorous 
analysis of climatologic, financial, electric power consump- 
tion and/or inventory change data for example, identifying the 50 
trends is very important. Under some circumstances non- 
periodic terms (i.e., the trends) may overwhelm the result, 
e.g., in computing the correlation function and spectral analy- 
sis, for example. In these instances it is equally important to 
remove the trend or detrend the data before arriving at mean- 55 
ingfi.il spectral results. Presently, because there is no precise 
mathematical definition for the aperiodic trend, trends are 
determined on a totally ad hoc basis. 

Typically, the trend is selected as the result of a moving 
mean, a regression analysis, a filtered operation or simple 60 
curve fitting with an a priori functional form. These trend 
approximations are all subjectively determined based on cer- 
tain idealized assumptions. Furthermore, normally, the data is 
not detrended using the same trend approximation. Instead, a 
typical trend determining method using a simple linear func- 65 
tion or straight line base is selected to define a detrend zero 
reference and the data is rezeroed by removing the reference. 
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Since especially for non-linear and non stationary data, the 
source of the particular trend is the very same mechanisms 
that generated the data, neither the trend approximation nor 
detrended data is particularly accurate. Consequently, 
approximating the trend using a linear fitting or with a moving 
mean makes little sense where the underlying mechanism is 
certainly nonlinear, non stationary and may lack a definitive 
time scale, e.g., for real world phenomena such as data from 
climatologic studies and financial data. 

Similarly, typical regression analysis and/or filtering are 
based on stationarity and linearity assumptions that may fit 
the data well, even for nonlinear regressions. However on a 
more fundamental level, a fortuitous regression analysis and/ 
or filtering fit for particular data, does not justify applying a 
particular regression formula globally as a time independent 
function. Thus, these various prior art curve fitting 
approaches use an a priori determined functional form that 
does not necessarily match the trend to the same underlying 
mechanisms, i.e., those that are inherent in the data. For 
example, financial data analysis pioneers, R. F. Engle and C. 
W. J. Granger produced market prediction models that are 
useful for a special class of non stationary processes. Engle 
and Granger treated the financial market as a special Auto- 
Regressive-Integrated-Moving-Average (ARIMA) process 
that is controlled by a series shocks and relaxations. Of 
course, as Engle and Granger have acknowledged, not all non 
stationary data satisfies their special assumptions. So, regard- 
less of the fit, it is very likely that the selected trend has a 
simplistic functional form that may not support the underly- 
ing mechanisms. Unfortunately, suitable analysis methods 
are not available for the vast majority of data from non sta- 
tionary and nonlinear signal. 

Thus, there is a need for a general method for identifying 
and determining trends from non stationary data and further, 
for detrending and determining the variability of such non 
stationary data. There is especially a need for a way to deter- 
mine the trend and to detrend data from non stationary and 
nonlinear processes that do not rely on extrinsic functional or 
simplifying assumptions. 

SUMMARY OF THE INVENTION 

It is an aspect of the invention to define the trend for non 
stationary and nonlinear signals; 

It is another aspect of the invention to intuitively and 
directly determine a precisely defined trend in any data 
including non stationary time varying data; 

It is yet another aspect of the invention to extract the trend 
from non stationary and nonlinear data signals such that the 
trend is defined by the same mechanisms from which the data 
are collected; 

It is yet another aspect of the invention to extract a trend 
that is an intrinsically fitted monotonic function within a data 
span, or a function in which there can be at most one extre- 
mum; 

It is yet another aspect of the invention to detrend non 
stationary time varying data simply by removing an intrinsi- 
cally determined trend and to determine the variability of the 
data after trend removal. 

The present invention relates to an apparatus, computer 
program product and method of analyzing non stationary 
time varying phenomena. A representation of a non stationary 
time varying phenomenon is recursively sifted using Empiri- 
cal Mode Decomposition (EMD) to extract intrinsic mode 
functions (IMFs). The representation is filtered to extract 
intrinsic trends, e.g., by combining a number of IMFs. The 
intrinsic trend is inherent in the data and identifies an IMF 
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indicating the variability of the phenomena. The trend also 
may be used to detrend the data. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The foregoing and other objects, aspects and advantages 
will be better understood from the following detailed descrip- 
tion of a preferred embodiment of the invention with refer- 
ence to the drawings, in which: 

FIG. 1 shows a block diagram example of a computer 
system including a computer suitable for programming with 
the inventive method; 

FIG. 2A shows a flow diagram example of iteratively 
extracting IMFs for any collection of data (X(t)); 

FIG. 2B shows a flow diagram example of extracting each 
IMF; 

FIG. 2C shows an example of Hilbert Spectral Analysis to 
derive trends and variability from the data and, also optionally 
detrend the data; 

FIG. 3 shows the NASDAQ daily closing index value sig- 
nal from its inception on 1 1 Oct. 1984 to 29 Dec. 2000; 

FIGS. 4A-I show the IMF components extracted from the 
from the NASDAQ data of FIG. 3 using CE(1000,5); 

FIGS. 5A-H show recombining the IMF components to 
reconstitute the original data signal; 

FIGS. 6A-B show a comparison of the intrinsically derived 
trend for the NASDAQ data with trends derived according to 
prior art methods and the same data detrended using those 
prior art trends; 

FIGS. 7A-D show identified temporally local time periods 
for partially reconstructed signals; 

FIGS. 8A-B show a comparison of short term trends deter- 
mined over a single year, 2000, using both traditional 30 and 
60 day moving averages and the analogous extracted local 
time period; 

FIGS. 9A-C show the time scale distribution for the 
extracted local time period of FIG. 7D and compare the 
detrending using 30 and 60 day moving averages with mov- 
ing averages determined from extracted local time period; 

FIGS. 10A-C show a comparison of NASQAD NASDAQ 
data variability as determined by prior art methods against 
variability according to a preferred an embodiment of the 
present invention. 

FIGS. 11A-B show the time scale detennination for the 
variability; 

FIGS. 12A-C show scattering plots of comparing the vari- 
ability measurements of FIG. 10A-C; 

FIG. 13 shows an example of a comparison of significance 
test results for the gain variability measurements and HL 
variability measurements with preferred EMD variability 
measurements; 

FIG. 14 shows an example of the results of summing the 
last six IMFs and their normalized values; 

FIGS. 15A-C show scattering plots of the relationship of 
the trends of FIG. 14 in more detail; 

FIGS. 16A-C show the variability values as a function of 
the market price slope or gradient of the statistically signifi- 
cant components and in particular, for examining higher vari- 
ability during the periods of downward market trend; 

FIG. 17 shows the annual global surface temperature 
anomalies from 1850-2003, i.e., the residual after removing a 
climatologic mean based on thirty years mean from 1961 to 
1990; 

FIGS. 18A-L show mean IMFs and the standard deviation 
of IMFs from each of ten different sets of sifting criteria, for 
the S from 4 to 13; 
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FIG. 19 shows the statistical significance test and indicates 
the significance of the last two IMFs; 

FIG. 20 shows the mean overall trend with its standard 
deviation; 

5 FIGS. 21A-B show selecting a trend with slightly shorter 
time scale than the overall trend by including the next com- 
ponent and the time scale for the trend; 

FIGS. 22A-B show selecting a trend with an even shorter 
time scale than the trend of FIG. 21A by including the next 
10 component and the time scale for the trend; 

FIGS. 23A-B show selecting a trend with an even shorter 
time scale than the trend of FIG. 22A by including the next 
component and the time scale for the trend. 

15 DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS 

Computer for Implementing Preferred Embodiment Meth- 
ods 

20 Turning now to the drawings and, more particularly, FIG. 1 

shows a block diagram example of a computer system 50 
including a computer 52 suitable for programming with the 
inventive method. The computer system 50 may include a 
keyboard 54 and mouse 56 for human interaction. When 
25 programmed with the inventive method, the computer 52 is 
analogous to an electrical filter or a mechanical sieve: it 
separates digital non stationary time varying signal data into 
series of Intrinsic Mode Functions (IMFs) according to their 
time scales. The digital data are filtered in a manner analogous 
30 to electrically filtering harmonics or a mechanical sieve 
which separates aggregated sand particles according to their 
physical size. An IMF represents an oscillation mode embed- 
ded in the data, as defined by the zero-crossings and involves 
only one mode of oscillation. 

35 Each IMF satisfies the following two conditions: (a) the 
number of extrema and the number of zero-crossings must 
either be equal or differ at most by one within the data, and (b) 
at any point, the mean value of upper envelope defined by the 
maxima and the lower envelope defined by the minima is 
40 zero. In other words, each IMF represents only one group of 
oscillation modes or time scales and no riding waves are 
allowed in an IMF. The first IMF condition is somewhat 
similar to the traditional narrow band requirements for a 
stationary Gaussian processes. Conceptually, the second con- 
45 dition modifies the classical global requirement to a local one. 
Furthermore, the second condition has the desirable result 
that the instantaneous frequency will not have unwanted fluc- 
tuations induced by asymmetric wave forms. Mathematically 
ideally, the second condition should be “the local mean of the 
50 data being zero.” For nonstationary data, a local time scale is 
necessary to compute the mean, which is otherwise impos- 
sible to define. Fortunately, the present invention obviates the 
need for defining a local time scale to fulfill the second con- 
dition. When applied to non stationary time varying signals to 
55 extract IMFs (e.g., from financial data representative of the 
daily NASDAQ index or geophysical data representative of 
annual global surface air temperature anomalies) trends are 
automatically identified and may further be removed to 
detrend the data. 

60 Thus, for application to non stationary time varying data, 
the local mean of the signal envelopes are used to force the 
local symmetry. The signal envelopes are defined by the local 
maxima and the local minima. This approximation avoids 
defining a local averaging time scale. With the preferred 
65 physical approach and approximation, an instantaneous fre- 
quency is extracted under most conditions. Moreover, even 
under the worst case conditions, the instantaneous frequency 
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so defined is still consistent with the nature of the system or 
phenomenon being studied and represents the system or phe- 
nomenon much more accurately than prior techniques such as 
those based on Fourier analysis. 

To facilitate analyzing non stationary time varying signals, 
the computer system 50 also includes an input device 58, 
and/or an imaging device 60, e.g., for collecting image infor- 
mation on a physical phenomenon and generating non sta- 
tionary time varying data representative thereof. Input device 
58 may include an antenna 61 connected to a receiver 63 for 
receiving transmitted data, for example, from the trading 
floor, or a remote weather station for climatological studies. 
The computer 52 may include an intrinsic mode extractor 5 1 , 
a low pass filter 53, a proto-IMF extractor 55 an IMF sifter 57 
and memory 70. Similarly, the computer system 50 also may 
include a display 62 or other output device such as a cathode 
ray tube or flat panel display, a printer 64 and any other 
suitable output device 66. Each of these outputs (62, 64. 66) 
maybe capable of generating or otherwise handling color 
outputs because, for example, the Hilbert Spectrum may be in 
color. The generalized output device 66 may also be con- 
nected through a network to the computer 52, e.g., in a local 
area network (LAN) or a wide area network (WAN). In this 
way, the non stationary time varying signals may be inputted 
from/transferred to the network. Thus, all outputs can be sent 
between remotely connected locations over such a network 
connection. 

Additionally, the computer system 50 also may include a 
mass storage device 68. The mass storage device 68 may be a 
hard disk, floppy disc, optical disc, etc. The mass storage 
device 68 may be used to store a computer program that 
performs the preferred method when loaded into the com- 
puter 52. More particularly, a computer program embodiment 
of the invention may be loaded from the mass storage device 
68 into the internal memory 70 of the computer 52. Loading 
the computer program embodiment transforms the general 
purpose computer 52 into a special purpose machine that 
implements the invention. Even more particularly, each step 
of a preferred embodiment method may transform at least a 
portion of the general purpose computer 52 into a special 
purpose computer module implementing that step. For 
example, when a sifting step is implemented on the computer 
52, the result is a computer implemented sifting apparatus 
(sifter) that performs the sifting functions of sifting step. 
Alternatively to storing the computer program embodiment 
locally, the input device 58 may be connected through a 
network connection to off-line storage that supplies the com- 
puter program to the computer 52. 

Other preferred embodiments of the invention include 
firmware embodiments and hardware embodiments wherein 
the inventive method is programmed into firmware (such as 
EPROM, PROM, flash memory or PLA) or wholly con- 
structed with hardware components. Constructing such firm- 
ware and hardware embodiments of the invention would be a 
routine matter to one of ordinary skill using known tech- 
niques. 

Article of Manufacture 

Still further, the invention disclosed herein may take the 
form of an article of manufacture, e.g., disc 72 in FIG. 1. More 
specifically, the article of manufacture may be any computer- 
usable medium, including a computer-readable program code 
embodied therein wherein the computer-readable code 
causes computer 52 to execute the inventive method. 
Examples of such a computer-usable medium 72 include a 
CD-ROM, a DVD-ROM, high density moveable storage, 
flash memory storage or a computer diskette. When the con- 
tents of the disc 72 are loaded into the mass storage device 68, 
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the computer-readable program code stored therein is trans- 
ferred into the computer 52. In this way, the computer 52 may 
be instructed to perform the inventive methods disclosed 
herein. 

5 

METHOD FOR CARRYING OUT THE 
INVENTION 

Rather than an extrinsically determined parameter, as 
io defined herein, the trend is an intrinsic data property and part 
of the signal generating the data and is the residue determined 
by this Empirical Mode Decomposition (EMD). For a suit- 
able example of EMD to determine IMFs and residue, see 
U.S. Pat. No. 5,983,162 entitled “Computer Implemented 
15 Empirical Mode Decomposition Method, Apparatus, and 
Article of Manufacture” issued Nov. 9, 1999 to the assignee 
of the present invention and incorporated herein by reference. 
See also, U.S. Pat. No. 6,782,124 entitled “Three Dimen- 
sional Empirical Mode Decomposition Analysis Apparatus 
20 Method and Article of Manufacture” issued Aug. 24, 2004 to 
Per Gloersen and assigned to the assignee of the present 
invention and incorporated herein by reference. The same 
mechanisms that generate the observed signal or measured 
data drive the intrinsic trend. More specifically, for an 
25 observed signal a trend exists only within a given data span 
(i.e., is temporally local) and is an intrinsically fitted mono- 
tonic function for that given data span (within the time of the 
data span), or a function in which there can be at most one 
extremum within the span. Accordingly, the trend is tempo- 
30 rally local and defined with reference to local time ranges, i.e., 
a local length, scaled in terms of data length. Further, the trend 
defined according to a preferred embodiment of the present 
invention has a more general meaning than the overall resi- 
due, because various local scales can be selected from within 
35 the data that clearly define a local trend. 

A definite local time scale must be determined in financial 
data, for example, to meaningfully model most financial sys- 
tems. Similarly, although it is important for societal and eco- 
nomical reasons to identify climatological trends, for 
40 example, such data typically is fraught with inherent time 
scaling variability. However, once intrinsic local periods or 
scales are defined, local trends may be extracted and the data 
may be subjected to meaningful analysis within the defined 
local scales. Thus, according to a preferred embodiment of 
45 the present invention, trends may be defined for any time 
series signal with varying local time scales, as long as a time 
scale may be defined within which the trend exists logically. 
Once the trend is defined, the data may be detrended simply 
by removing the trend and the variability can be extracted. 

50 Generally, as canbe seen from the flow diagram example of 
FIG. 2A, a non stationary time varying signal that may be 
non-linear is subjected to the Hilbert-Huang Transform 100 
to extract intrinsic trends and variability data from the signal 
and, optionally detrend the signal. So, for any collection of 
55 data 102 representing the non stationary time varying signal 
(X(t)), IMFs are iteratively extracted 104, 106. Essentially, 
each IMF represents a single oscillation mode embedded in 
the data, defined by zero-crossings. In step 108, the EMD 
components are subjected to Hilbert Spectral analysis to 
60 adaptively extract and derive trends and variability from the 
data and, also detrend the data according to a preferred 
embodiment of the present invention. This is in contrast to 
prior art approaches, which typically define an extrinsic 
trend, selected for example, using pre-selected simple func- 
65 tional forms and time periods. Moreover, the intrinsic trend 
can be selectively locally valid within part of data, i.e., less 
than a full local wavelength. Where there were two or more 
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extrema of the same sense, the data may represent a cycle 
instead of a trend. Thus, the present invention avoids the well 
known confusion in that “one economist’s ‘trend’ can be 
another’s ‘cycle.’” 

FIG. 2B shows a flow diagram example of extracting IMFs 
in step 104 according to a preferred embodiment of the 
present invention. As indicated hereinabove, each IMF satis- 
fies the following two conditions: (a) the number of extrema 
and the number of zero-crossings must either be equal or to 
differ at most by one within the data, and (b) at any point, the 
mean value of the upper envelope defined by the maxima and 
the lower envelope defined by the minima is zero. In other 
words, each IMF represents only one group of oscillation 
modes or time scales and no riding waves may be included in 
an IMF. It should be noted that each of these local time scales 
are the bounds of a temporally local trend. 

So, beginning in step 1042 local extrema are identified in 
the signal X(t) or the partially decomposed signal X,(t), i.e., 20 
after at least one has been removed. Then, in step 1044 the 
signal envelope is marked at the signal extrema to identify 
local maxima which defines an upper envelope edge and local 
minima which defines a lower envelope edge. Each edge may 
be identified independent of or simultaneously with the other. 25 
In step 1046 the mean (designated as m,) of the envelope is 
determined as the mean of the upper and lower edges. In step 
1048 as the difference between the data or partially sifted data 
and the mean from step 1046 determines a proto-IMF (h,), 30 
i.e., X,(t)-m,=h,. Since with each IMF extraction, previously 
hidden data characteristics may become more pronounced 
and so, become a local extremum, the proto-IMF is checked 
in step 1050 to determine if additional sifting has become 
apparent after each extraction. If additional data characteris- 35 
tics must be removed, then the proto-IMF is treated as the data 
in step 1052 and local extrema are identified in step 1042. 
Ideally, however, the proto-IMF, h,, is the IMF in step 1054. 
Thus, in step 104 of FIG. 2A, sifting is repeated until no IMFs 
remain to be extracted, i.e., the remaining data contain only 40 
one extremum and is the residue defining the overall trend. 

It is this sifting of the proto-IMFs that eliminates riding 
waves and is necessary for separating the final IMFs 1054 and 
for defining a meaningful instantaneous frequency. Further, 4 _ 
since the riding waves are eliminated, the wave profiles are 
symmetric, which may be necessary if the disparity between 
neighboring wave amplitudes is too large. So, the proto-IMF 
is repeatedly sifted as many times as necessary to eliminate 
any remaining riding waves and improve IMF symmetry. This 50 
repeated sifting may be used to reduce the extracted signal 
from a proto-IMF to an IMF by treating each proto-IMF as the 
data, and then, applying 

K-”>, i=h„ 55 

used repeatedly sifting up to k times to yield 

hi(k-T)- m ik = Kk> 

where h ik is the i‘ h IMF component designated as c, h, 4 . Over- 60 
all, c, should contain the finest scale or the shortest period 
component of the sifted signal and the first component, for 
example, can be separated from the rest of the data by 

X{t)-c l =r l , and generally thereafter, 65 


r^f-c^r^ and in summary for n components, 


X(t) = Y J c i +r, 


Each of these empirical mode components usually is physi- 
cally meaningful within characteristic scales that are self- 
defined by the physical data and the residue r, defines the trend 
within those characteristic scales. As long as the residue, r„ 
contains periodic components in step 106 of FIG. 2A, it is 
treated as the data and sifted in step 104 as described above 
with reference to FIG. 2B until all of the components have 
been removed. The overall trend r n is the last extracted com- 
ponent and is identified when sifting produces a constant, a 
monotonic function or a function with only one maximum 
and one minimum front which additional IMFs camiot be 
extracted. Since even if the data X(t) has a zero mean, its final 
residue may be non zero and the final residue defines the 
overall trend. Thus, the signal X(t) decomposes into n-ern- 
pirical mode components including IMF components c, and 
the final residue r„. These EMD components (c,) are subjected 
to Hilbert Spectral Analysis in step 108 to derive trends and 
variability from the data and, also optionally detrend the data 
according to a preferred embodiment of the present invention. 

FIG. 2C shows an example of the Hilbert Spectral Analysis 
step 108 to derive trends and variability from the data and, 
also optionally detrend the data. First in step 1080, the Hilbert 
Transform is applied to each component (c,.). Then, in step 
1082, the data are expressed as a real part. In step 1084, the 
Hilbert Spectrum is extracted from the real part . In step 1 086, 
the marginal spectrum is extracted front the Hilbert Spec- 
trum. In step 1088 variability and trend filtering frequencies 
are identified from the spectrum data. In step 1 090, the signal 
is low passed filtered recombining a number of lower fre- 
quency components with the final residue r„ to generate inter- 
mediate trends 1 092 which are either a constant or monotonic 
function. Optionally in step 1094, the trend 1092 may be 
subtracted from the signal to detrend the data. Similarly the 
cut off frequency identified in step 1088 identifies corre- 
sponding IMF components that indicate the variability of the 
signal at those trends. So, optionally in step 1096 the corre- 
sponding IMF components may be rectified by taking the 
absolute value of the IMF and, normalized to the original 
signal. Thus, the corresponding IMF components, optionally 
rectified and normalized, are the corresponding variabilities 
1098. 

In step 1080 the Hilbert transform, Y(t), is computed and 
applied to each IMF component to compute the instantaneous 
frequency in, i.e., applying Hilbert spectral analysis to EMD. 
In particular, the Hilbert transform, which is the convolution 
of X(t) and 1/t, has the form: 


Y(t) = - 


A pf* 

n J t 


X(t') 


dr'. 


where P is the Cauchy principal value. The Hilbert transform 
emphasizes the local component of X(t), even though the 
transform is global . It is well known that the Hilbert transform 
exists for all functions of L p class and so, X(t) and Y(t) form 
a complex conjugate pair of an analytic signal, Z(t). In polar 
coordinates: 
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Z(t)=X(t)+iY(f)=a{ t)e m< -‘\ where 
a (t)=[X 2 (t) + iY 2 (t)] 1 ' 2 ; 

0 (?) =arctan ( Y{t)!X(t)) . 

The polar coordinate expression fiirther clarifies the local 
nature of this representation, providing a local fit of an ampli- 
tude and phase varying trigonometric function to X(t). Using 
the Hilbert transform, the instantaneous frequency may be 
defined as: 

(0 (t)=dQ(t)ldt. 

So, once the Hilbert transform is obtained for each IMF, in 
step 1082 the original data can be expressed as the real part, 
RP, having the form: 


X '(!) = RP^ a j (t)e , to wd ' . 
j= 1 


The Hilbert transform yields both an amplitude and fre- 
quency of each component as functions of time and has an 
energy that is the component’s Hilbert Amplitude Spectrum, 
H(w,t), in time-energy-frequency space or, simply, its Hilbert 
Spectrum in 1084. Although the final residue (the trend, r„) 
has been omitted, if physical considerations justify inclusion, 
the Hilbert transform can treat the monotonic or constant 
trend as part of a longer oscillation. However, normally, it is 
omitted to avoid the possibility that energy in the residual 
trend would overpower the rest of the spectrum. In particular 
any benefit derived from including the trend component sel- 
dom outweighs the uncertainty introduced by any increase in 
the time period to include the trend. Further, normally, inter- 
est lies in oscillatory components including the information 
contained in the other low energy but clearly oscillatory com- 
ponents. 

It should be noted that expanding the same signal in a 
Fourier series representation has the form: 


h{oj) = I H(co,t)dr. 
Jo 


Generally, the marginal spectrum represents the cumulated 
amplitude over the entire data span in a probabilistic sense. 
From this representation, the areas of interest, e.g., changing 
trends and periods of intense variability, may be identified in 
step 1088. 

It should be noted that, in particular, traditional signal 
filtering may be applied to the IMF components, especially 
for either nonlinear or nonstationary signals or both. Advan- 
tageously, filtering the IMF components from such signals 
that are both nonlinear and nonstationary avoids wave form 
deformation and distortion from inadvertently generated har- 
monics of all ranges including those outside of the filtering 
range and that prior art approaches have frequently encoun- 
tered. 

Accordingly where periods T of interest identified by the 
Hilbert spectrum peaks of a filtering point as the k ,; ’, in step 
1090 and 1094 the signal is easily filtered with standard 
time-space filtering, e.g., low pass filtering k of the n IMF 
components of the signal to extract a trend can be simply 
expressed as 


Xu (0 = 2 cj + r„-, 

t. 


high pass filtering the first expressed as 


k 

XhkU) — Cj\ 


10 

contribution can be measured for each frequency in step 1086 
using the marginal spectrum h(co), which has the form 


5 r T 


X(t) - RP jr Ojc'H'" 
M 


and 

4; ' band pass filtering expressed as 


with both and to, term being constants rather than time 
dependent functions. So, the IMFs are a generalized Fourier 
expansion of the signal. In contrast to the constant value 
Fourier representation, the preferred embodiment IMF 
expansion including the above real part representation exhib- 
its greatly improved efficiency with variable amplitude and 
instantaneous frequency components that accommodate non- 
linear and nonstationary data. Further, the IMF expansion 
clearly separates amplitude from frequency modulations, 
overcoming the Fourier expansion restriction of selecting 
constant amplitudes and fixed frequencies to arrive at a vari- 
able amplitude and frequency representation in step 1084, 
i.e., the Hilbert Spectrum. 

Further, to represent energy density, in step 1086 squared 
values of amplitude may be substituted to produce the Hilbert 
Energy Spectrum as well. However, for more quantitative 
results the skeleton Hilbert Spectrum may be presented. 
Additionally, two-dimensional smoothing (e.g., the Laplace 
Transform) may be applied to the skeleton for a “fuzzy” or 
“smeared” view. In particular, the total amplitude (or energy) 


k 

XbkU) ~ ^ \ Cj. 
50 b 


Advantageously, time-space filtering according to a preferred 
embodiment of the present invention preserves the lull non- 
55 linearity and nonstationarity in physical space of each filtered 
signal including both the trend 1092 and the high passed 
filtered result of step 1094. 

Optionally, the absolute value may taken of the high pass 
filtered result and normalized to the corresponding instanta- 
60 neons signal X(t) and variability V(t;T) may be selected to 
have the form: 


65 


V(r; T) = 


l-y«(r)| 
XU) ■ 
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Variability is somewhat different than the more traditional 
market analysis parameter of volatility. Instead, while vari- 
ability 1098 is a function of time, it also is time scale depen- 
dent, i.e., dependent upon T. Further, variability 1098 is a 
simple and direct measure of the market volatility and is a unit 
free parameter because it is a normalized to the market value. 

Once the trend 1092 is generated, detrending is a relatively 
simple operation. The data X(t) is detrended by subtracting 
the trend, i.e., X(t)-X a (t). 

Advantageously, intrinsic trends may be extracted for any 
non stationary time varying signal according to a preferred 
embodiment of the present invention and further, the relation- 
ship of the trend with the IMFs may be described. Moreover, 
the preferred embodiment trend identification approach also 
may be applied to determining various trends with nuanced 
time scales, including to real world non stationary time vary- 
ing signals, such as financial data or climatological data. 
Thus, a better understanding and appreciation of the advan- 
tages and uses of the present invention may be achieved 
through specific examples provided hereinbelow of applica- 
tion to both financial data (e.g., the NASDAQ daily index) and 
climatological data (e.g., the annual global surface tempera- 
ture data). It should be noted that in both specific examples, 
the data may not be representative of general principals, i.e., 
the daily NASDAQ index may not represent the dynamics of 
the financial market well; and, anomalies in the annual global 
surface temperature data may not portray climate change 
mechanisms clearly. However, data for these two real world 
phenomena are readily available and serve well to illustrate 
the present invention. The present invention has application to 
determining various trends with nuanced time scales in any 
type of data. 

The Daily NASDAQ Index 

So, for example, FIG. 3 shows the NASDAQ daily closing 
index value signal 110 from its inception on 1 1 Oct. 1984 to 
29 Dec. 2000, a total of 4099 days. FIGS. 4A-I show the IMF 
components 112, 114, 116, 118, 120, 122, 124, 126, 128, 
respectively, extracted as described in FIG. 2B from the from 
the NASDAQ data of FIG. 3 using CE(1000,5) (extrema 
sifting with stoppage S number equal to 5 consecutive proto- 
IMFs, i.e., with same and equal number of extrema and zero- 
crossings). So, the overall trend in this example is the 9 th 
residue 128 of FIG. 41. 

FIGS. 5A-H show step-by-step 130, 132, 134, 136, 138, 
140, reconstituting the original data from the IMF compo- 
nents 112, 114, 116, 118, 120, 122, 124, 126, 128 of FIGS. 
4A-I, with reference to the original signal 110. So, FIG. 5A 
shows the original signal 110 overlaying the trend 128. As 
each subsequent IMF component 126, 124, 122, 120, 118, 
116, 114 is added, the reconstituted signal 130, 132, 134, 138, 
140 approaches the original until the first IMF 1 12 is added to 
completely reconstitute the original 110 in FIG. 4H. Of 
course, the order in which the selected IMFs are combined 
112, 114, 116, 118, 120, 122, 124, 126, 128 is unimportant as 
the result will be the same, i.e., the original signal 110. 

FIGS. 6A-B show a comparison of the intrinsically derived 
overall trend 128 for the NASDAQ data 110, derived accord- 
ing to a preferred embodiment of the present invention, with 
overall trends 150, 152 derived according to prior art meth- 
ods; and detrended 154, 156 using those prior art overall 
trends 150, 152 and detrended 158 according to a preferred 
embodiment of the present invention. Thus, extrinsically 
determined trends 150, 152, which are determined linearly 
and exponentially, respectively, are compared with the final 
residue 128, which is the last IMF component of FIG. 41, 
overlaying the original NASDAQ data signal 110. By com- 
parison the intrinsically determined overall trend 128 pro- 
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vides much more valuable information and more closely pre- 
dicts the signal for the whole data span. Further, comparing 
the detrended results of FIG. 6B shows that the traditional 
detrended data 154, 156 are much poorer results than pre- 
5 ferred embodiment detrended data 158. It should be noted 
that typical current stock market models use an exponential fit 
152, 156 which, although provides a much better fit than the 
linear fit 150, 154, still falls far short of the preferred embodi- 
ment EMD fit 128, 158. Similarly, the variability statistics of 
to the fittings is much improved for a preferred embodiment 
detrending than from prior art methods. Specifically, it may 
be shown that the mean of the variability is zero for a linear fit, 
73.1187 for an exponential fit and 0.3588 for a preferred 
EMD fit; and, the standard deviation of the variability is 
15 559.09 fora linear fit, 426.66 for an exponential fit and 238. 10 
for a preferred EMD fit. 

Further, it may be shown that using traditional polynomial 
curve fitting techniques, the NASDAQ data 110 may be mod- 
eled with a 4 ,; ' power polynomial. For one such polynomial 
20 model, the mean deviation is only -0.0079 and the standard 
deviation is 255.36, which is only slightly larger than 238. 10 
for the preferred EMD fit. It should be noted, however, that 
other than that this A ,h power polynomial is a better fit than 
either of the simple linear and exponential fittings, there is no 
25 other justification for its application to market, because there 
is nothing to indicate that such market data is time dependent 
or derived from time dependent processes that are fourth 
order functions of time. 

Generally, the trend has meaning for the entire time period 
30 of interest, but is not particularly helpful in extracting other 
useful data. Short range trends overmuch finer time scales are 
necessary to determine meaningful temporally local param- 
eters and behaviors, e.g., a weighted mean and a standard 
deviation. So, while the final residue 128 provides an overall 
35 trend for the total period for the signal 110, that may not be 
particularly useful for understanding and determining these 
local or short range trends. Short range trends within shorter 
time periods may be determined by identifying an upper 
bound of by the variability zero-crossings within the range of 
40 interest. In particular, the range or time scale may be deter- 
mined by the time spans between a various combinations of 
the critical points defined as the union of all the zero-cross- 
ings and extrema. 

So, for a particular time range, the most local time scale and 
45 also the finest resolution, is the distance between the neigh- 
boring critical points, i.e., between an extrema to the next 
zero-crossing. Next within the same time range is a half wave 
cycle, i.e., the time either between two consecutive extrema (a 
minimum to the next maximum), or between two consecutive 
50 zero-crossings (up zero-crossing to a down-zero-crossing). 
The full wavelength is the longest time scale within the time 
range and is also the most physical time scale. The full wave- 
length may be measured from one maximum (minimum) to 
the next maximum (minimum), or one up (down)-zero-cross- 
55 ing to the next up (down)-zero-crossing. Thus, for any given 
time range, a number of local time scales may be determined 
for representing different temporally local properties, i.e., 
within that given time range. 

FIGS. 7A-D show the local periods for partially recon- 
60 stnicted or low pass filtered signals. Combining the last two 
IMF components (126, 128 in FIGS. 3H and I as shown in 
FIG. 5B) improved data fit over the overall trend (FIG. 5A) 
because as can be seen in FIG. 7A, the filtered signal has an 
associated time varying period 160 that is much finer than the 
65 overall trend 128. Instead, a local trend may be identified in 
each of these periods. Further, the time varying period 160 
varies between 400 to 1200 days with a mean at near 800 
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days, based on the local data variability rates, which is an 
indication of the dynamics of the market. Accordingly, the 
mean value, which is 796.6 days or nearly 4 years (there are 
only 200 trading days in each year), might be considered 
meaningless for this data. However, using the mean value as 5 
a guide for the overall time scale, this four-year cycle is 
coincides with the time scale of the presidential election. This 
same four-year cycle has also appeared in the 30 years mort- 
gage interest rate studies. Thus, as might otherwise be 
expected, it appears possible/likely from FIG. 7A that politics 10 
influences the financial market. 

Similarly, in FIG. 7B, local trends can be defined for finer 
periods by opening the low pass filter slightly, including the 
last three IMF components 124, 126, 128 for an even better 1 . 
data fit of the as shown in FIG. 5C. In this example, the local 
period 162 varies between 200 to 600 days, with a mean of 
425.7 days or slightly over two years. Likewise, the local 
period 164 of FIG. 7C with the low pass filter including the 
last four IMF components 122, 124, 126, 128 varies between , f) 
80 and 420 with a mean value of 1 96; and, in FIG. 7D the local 
period 166 from including the last six IMF components 118, 
120, 122, 124, 126, 128 varies between 8 and 77 with a mean 
of 35.6 days. It should be noticed that the time scales are 
variable and process dependent and, in particular not deter- 
mined subjectively. Therefore, time scale selection may be 
effected only by selecting different IMF combinations, i.e., 
low pass filter selection. However, unlike prior art 
approaches, once the IMFs are selected, the data determine 
the time scale and, accordingly, the trend. Also the periods are 3(j 
determined by which successive IMF component are 
included and reveal nearly factor 2 change, indicating that the 
EMD filtered components may serve as a dyadic filter. 

FIGS. 8A-B show a comparison of short term trends deter- 
mined over a single year, 2000, using traditional 30 and 60 35 
day moving averages 1 70, 172, respectively, and based on the 
analogous extracted local time period 166 of FIG. 7D as 
determined according to a preferred embodiment of the 
present invention. For nonstationary financial data with no 
discemable local time scale and/or nonlinear processes and in 40 
the financial market analysis in particular, the most popular 
ways of defining trends has been regression analysis, moving 
averages or means. However, frequently stationary and, 
sometimes even, linear assumptions that are used for regres- 
sion analysis and filtering have proven problematic. Simi- 45 
larly, a pre-determined time length scale, e.g., 30 days and 60 
days in this example, must be selected to determine the mov- 
ing mean. Consequently, both of these powerful methods 
have proven unsatisfactory for the nonlinear and nonstation- 
ary processes. By contrast, however, as can be seen by the 50 
comparisons below, trends extracted according to a preferred 
embodiment of the present invention are truly adaptive. 

So, as can be seen from the example of FIG. 8A, it is 
difficult to justify selecting the 30 and 60 day fixed period or 
any fixed period for application to the whole non stationary 55 
time varying signal. Further, although adaptive, selecting the 
appropriate time scale is problematic in determining the mov- 
ing average. Even if an appropriate time scale were selected, 
it is impossible to know the future value. So, the moving mean 
170, 172 is derived from past data and is the backward mean, 60 
out of phase with the data. By contrast, the preferred embodi- 
ment EMD approach determines the intrinsic trend in real 
time and so, extracts the appropriate time scale from the data. 
FIG. 8B shows an example wherein the 30 day mean 170' and 
60 mean 172' are shifted by half of the averaging time dura- 65 
tion (i.e., 1 5 days, 30 days) and compared against a preferred 
embodiment intrinsic trend 174, which is in real time and 
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unshifted. Thus, the intrinsically fit trend 174 extracted 
according to a preferred embodiment of the present invention 
more closely tracks the data. 

FIGS. 9A-C show the time scale distribution 180 for the 
extracted local time period 166 of FIG. 7D and compare the 
detrending using 30 and 60 day moving averages with mov- 
ing averages determined from extracted local time period 
166. Thus using an adaptive time scale, the preferred embodi- 
ment overall fit to the data is better as shown by the detrended 
30 day moving average signal 181, 60 day moving average 
signal 182 and EMD detrended signal 183 of FIG. 9B and 
their distribution given in FIG. 9C. In particular, the Standard 
Deviations for the 60-day mean 184 is at 54.05; the 30-day 
mean 186 is at 40.38, and the preferred EMD mean 188 is at 
35.86 even though the mean period in this fitting is slightly 
longer than 30 days. Thus, the preferred embodiment data 
fitting provides a much truer fit. 

As noted hereinabove, the variability is a measure of the 
volatility of the financial data. It is understood that as used 
hereinbelow and for simplicity of di scussion, variability is the 
same as volatility for financial data, although variability is 
more appropriate and meaningful generically for other types 
of data such as climate change. Further, with regard to finan- 
cial data a traditional definition of variability is a measure of 
the market movement or the day-to-day market fluctuation. A 
common such measure of market movement is gain, G„ 
which is a discrete logarithm at the time t of the ratio of 
consecutive index values, S„ and S,.!, orthe logarithmic ratio 
of consecutive values (LCRV) and is defined as 


s, 

G, = log-2- . 
*->?—! 


Thus, 

G ,= log 5,-log 5 m . 

FIGS. 10A-C show a comparison of NASDAQ data vari- 
ability as determined by prior art methods against variability 
determined according to a preferred embodiment of the 
present invention, e.g., the residue after removing the trend. 
In particular, FIG. 10A shows an example of variability 
extracted as the gain of the NASDAQ signal 110. FIG. 10B 
shows an example of variability extracted as the difference 
between the highest and lowest values (HL) of the price of the 
NASDAQ signal 110 within one day. By contrast, FIG. IOC 
shows a measure of the variability as the first IMF component 
112 in FIG. 4B from EMD. Thus, it can be seen that the 
preferred embodiment variability measurement provides an 
indication of the fastest or the shortest periods of price fluc- 
tuation. 

It should be noted that the three methods do not afford a 
strict apples to apples comparison because, for example, the 
HL values are always positive, while both the gain and EMD 
values can be both positive and negative. So, to facilitate 
comparison, for example, the gain and the EMD variability 
may be rectified by taking their absolute values. Alternately, 
since the gain measures a ratio, the HL and the EMD may be 
normalized by the local daily price for a better comparison as 
in steps 1094-1098 of FIG. 2C. Additionally, the time scales 
are not common for all three variability examples of FIGS. 
10A-C. Both the gain and HL measurements are day-to-day 
or daily variability measurements. By contrast the preferred 
embodiment EMD variability measurement is the extracted 



US 7,464,006 B1 


15 

fluctuation from an underlying mean of the next IMF com- 
ponent 114. Therefore, the time scale of this underlying mean 
must be identified. 

FIGS. 11A-B show the time scale determination for the 
variability according to a preferred embodiment of the 5 
present invention. In particular, FIG. 1 1 A shows period deter- 
mination for the first IMF component 112, determined analo- 
gously to FIGS. 7A-D. The time period for this first IMF 
component 112 has a mean of 8.38 days and a range of 2-25 
days. FIG. 11B shows the distribution of the time scale of the to 
time periods of FIG. 11A. So, the preferred embodiment 
EMD variability measures essentially the price fluctuation 
over an 8-day mean period. 

FIGS. 12A-C show scattering plots of comparing the vari- 
ability measurements of FIGS. 10A-C. FIG. 12A shows a 15 
pairwise comparison of the gain and EMD variability mea- 
surements. FIG. 12B shows a pairwise comparison of the gain 
and HL variability measurements. FIG. 12C shows a pairwise 
comparison of the EMD and HL variability measurements. 
Both the EMD components and the HL variability measure- 20 
ments are derived from data extrema. So, as might be 
expected, these scattering plots indicate a lack of interrela- 
tionship except some hint of correlation between the HL and 
EMD measurements in FIG. 12C. Both reflect the high-low 
differences, although the EMD variability is over a period 25 
longer than one day, i.e. the major difference is the time scale. 
The fact that there is some correlation between the daily 
fluctuation and the fluctuation over a mean of 8 days suggests 
that the market has some persistence. 

FIG. 13 shows an example of a comparison of significance 30 
test results for the gain variability measurements 190 and HL 
variability measurements 192 with preferred EMD variability 
measurements 194. Clearly, all three 190, 192, 194 deviate 
from the white noise model. Thus, although it has been sug- 
gested that financial data are fractional Gaussian noise, the 35 
NASDAQ data of this example clearly does not comply with 
the fractal Gaussian test. 

FIG. 14 shows an example of the results of low pass filter- 
ing the last six IMFs 118, 120, 122, 124, 126 and 128 in FIGS. 
4D-I, and plotting their normalized variability values based 40 
on a line equal to three times the period as a statistical sig- 
nificance boundary indication. FIGS. 15A-C are scattering 
plots showing the relationship of the trends of FIG. 14 in more 
detail. FIG. 15A shows a pairwise comparison of the gain and 
EMD trends. FIG. 15B shows a pairwise comparison of the 45 
gain and HL trends. FIG. 15C shows a pairwise comparison 
of the EMD and HL trends. Thus, all methods show clear 
trend correlations. It should be noted that the trend values may 
be viewed as trends or variability values, depending on the 
point of view, i.e., whether the period of interest is shorter 50 
than the time scale adopted (a trend), or longer than the time 
scale (variability). 

Finally, FIGS. 16A-C show the variability values as a func- 
tion of the market price slope or gradient of the statistically 
significant components and in particular, higher variability 55 
during the periods of downward market trend may be exam- 
ined. So, FIG. 16A shows market price gradient for gain 
defined variability. FIG. 16B shows market price gradient for 
HL defined variability. FIG. 16C shows market price gradient 
for EMD defined variability. The tendency of having the 60 
higher variability values in the negative market gradient 
appears to be very clear. 

Global Surface Temperature Anomaly Data 

Global warming has become a contentious economic and 
political problem as well been as a scientific issue. Data on the 65 
annual global surface temperature anomalies have been 
painstakingly assembled and are available from the Climate 
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Research Unit (CRU), University of East Anglia, UK, at the 
Website, e.g., www.cru.uea.ac.uk/cm/data/temperature/. 
Although, there is no longer any doubt that there has been a 
warming trend over the last hundred years, the causes as well 
as the warming rates and, in particular, the precise trend 
remain unsettled. Studying the intermediate term variability 
of the climate may help to resolve some of these issues and, in 
particular, identify the climatological trends. 

FIG. 17 shows annual global surface temperature anoma- 
lies 200 from 1850-2003, i.e., the residual afteraclimatologic 
mean based on a thirty year mean from 1961 to 1990 had been 
removed. As defined according to a preferred embodiment of 
the present invention, the trend and the associated time scales 
are defined using EMD. After these data were analyzed using 
curvature sifting for the extrema sifting, too many inflection 
points remained, which suggested imbedded hidden scales. 
FIGS. 18A-L show mean IMFs 202, 204, 206, 208, 210, 212 
and standard deviation of IMFs 214, 216, 218, 220, 222, 224 
from each of ten different sets of sifting criteria, for the S 
number from 4 to 13. Each sifting gives six IMF components 
202, 204, 206, 208, 210, 212and214, 216, 218, 220, 222, 224, 
respectively, indicating that the siftings are producing 
extremely stable results. The standard deviation values 214, 
216, 218, 220, 222, 224 are about one order of magnitude 
smaller than the mean values 202, 204, 206, 208, 210, 212. 

FIG. 19 shows the statistical significance test results 230 as 
compared with white noise 232, which indicates that the last 
two IMFs 210, 212 are certainly significant. Accordingly, for 
this data it would not prove fruitful to pursue the trend to any 
finer scale, i.e., a time scale less than 20 years is hardly 
significant. FIG. 20 shows the mean overall trend 212 with its 
standard deviation 224. As with the above financial data, this 
trend is neither linear nor quadratic, but intrinsically deter- 
mined and as indicated by the 95% limit, highly reliable. FIG. 
21A shows a first low pass filtering with the next component 
210 included for trend 234 with a slightly shorter time scale 
and FIG. 21B shows the time scale 236 for the trend 234. The 
time scale 23 6 is certainly not constant, but it has a mean value 
slightly higher than 65 years. FIG. 22A shows an expanded 
low pass filter that includes the next component 208, for a 
trend 238 with yet a shorter time scale and FIG. 22B shows 
the time scale 240 for the trend 238. Here the mean value is 
dropped to about 20 years. FIG. 23 A shows a further 
expanded low pass filter that includes the next component 206 
for trend 242 with yet a shorter time scale and FIG. 23B shows 
the time scale 244 for the trend 242. More reliable trends are 
based on the longer time scale. Thus, other than the familiar 
overall global wanning trend, the 65 year cycle really stands 
out. 

Advantageously, a precisely defined trend may be intu- 
itively and directly extracted according to the present inven- 
tion from any data including non stationary time varying data. 
Once extracted, the trend serves for precisely detrending the 
data and determining data variance, both important data 
analysis steps. Instead of the various ad hoc extrinsic methods 
that have been used to determine the trend, or to detrend the 
signal source, the preferred embodiment method extracts an 
intrinsic trend from any nonlinear and nonstationary signal or 
time series. Thus, the extracted trend is determined by the 
same mechanisms that generate the data and, therefore, the 
trend is part of the data. In particular, the extracted trend is an 
intrinsically fitted monotonic function within the data span, or 
a function in which there can be at most one extremum. 
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Further, the preferred embodiment trend identification 
method adaptively extracts the intrinsic trend. Moreover, the 
extracted trend is more than simply the overall residue, but 
also defines an associated time scale. Once the trend is 
defined, the data may be detrended simply by removing the 5 
trend. The variability can be extracted as the residue of the 
data after trend removal and optionally normalized for com- 
parison against other variability model results. 

Furthermore, thepreferred signal decomposition method is 
applicable for extracting a trend from any data with different to 
simple intrinsic oscillation modes. The trend may be 
extracted from data that may have many different coexisting 
oscillation modes, at any given time superimposed on each 
other. Each mode, which may or may not be linear, has the 
same number of extrema and zero-crossings and oscillates 15 
symmetrically with respect to a local mean. 

While the invention has been described in terms of pre- 
ferred embodiments, those skilled in the art will recognize 
that the invention can be practiced with modification within 
the spirit and scope of the appended claims. 20 

What is claimed is: 

1. A computer implemented method of extracting trend 
data from a non- stationary time varying phenomena, said 
method comprising the steps of: 

inputting a representation of a non-stationary time varying 25 
phenomenon; 

recursively sifting said representation using Empirical 
Mode Decomposition (EMD) to extract an intrinsic 
mode function (IMF) indicative of an intrinsic oscilla- 
tory mode wherein said IM-F is a proto-IMF until said 30 
proto-IMF has a zero mean; 

filtering said recursively sifted representation including the 
step of identifying cut off frequencies; 
wherein identified said cut off frequencies identify one or 
more IMF indicating variability of said non-stationary' 35 
time varying phenomena; 
rectifying said identified one or more IMFs; 
normalizing rectified said identified ones; 
low pass filtering said representation responsive to said cut 
off frequencies; extracting trend data from said filtered 40 
representation; and 
displaying said extracted trend data. 

2. The computer implemented method of claim 1, wherein 

the step of low pass filtering comprises partially reconstruct- 
ing said representation. 45 

3. The computer implemented method of claim 2, wherein 
said partially reconstructed representation has the form 


X*(0 = Cj + r„, 


where k identifies a cut off frequency, n is the number of 
EMD components, c 1 represents IMFs, r„ is the final 55 
residue and X ft (t) is the extracted trend. 

4. The computer implemented method of claim 1, wherein 
the step of identifying said cut off frequencies further com- 
prises the steps of: 

applying a Hilbert transform to each said IMF ; 60 

expressing said representation as a real part; 
deriving the Hilbert spectrum from said real part; and 
deriving a marginal spectrum from said Hilbert spec- 
trum, energy in said marginal spectrum identifying 
said cut off frequencies. 65 

5. The computer implemented method of claim 1, wherein 
said representation is non-linear. 
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6. A computer program product for extracting trend data 
from a non-stationary time varying phenomena, said com- 
puter program product comprising a computer readable 
medium having computer readable program code compris- 
ing: 

computer readable program code for receiving a represen- 
tation of a non-stationary time varying phenomenon; 
computer readable program code for extracting intrinsic 
mode functions (IMFs) for Empirical Mode Decompo- 
sition indicative of an intrinsic oscillatory mode from a 
received said representation wherein at least one of said 
IMFs is a proto-IMF; 

the Empirical Mode Decomposition comprises computer 
readable program code for recursively sifting said rep- 
resentation comprising computer readable program 
code for extracting a proto-IMF and computer readable 
program code for extracting said IMF from said proto- 
IMF by recursively sifting said proto-IMF until said 
proto-IMF has a zero mean, said recursively sifted 
proto-IMF being said extracted IMF; 
computer readable program code for selectively filtering 
said extracted IMFs in said representation, said filtered 
representation indicating trends in said non-stationary 
time varying phenomenon; 

the computer readable program code for selectively filter- 
ing comprising: 

computer readable program code for identifying cut off 
frequencies, 

wherein identified said cut off frequencies identify one or 
more IMF indicating variability of said non-stationary 
time varying phenomena; 

the computer readable program code for identifying said 
cut off frequencies further comprising: 

a) computer readable program code for applying a Hil- 
bert transform to each said IMF; 

b) computer readable program code for expressing said 
representation as a real part; 

c) computer readable program code for deriving the 
Hilbert spectrum from said real part; 

d) computer readable program code for deriving a mar- 
ginal spectrum from said Hilbert spectrum, energy in 
said marginal spectrum identifying said cut off fre- 
quencies; 

e) computer readable program code for rectifying said 
identified one or more IMFs; and 

f) computer readable program code for normalizing rec- 
tified said identified ones; 

computer readable program code for low pass filtering said 
representation responsive to said cut off frequencies, 
said low pass filtered representation being said extracted 
trend; and 

an output device for displaying said filtered representation 
indicating trends. 

7. The computer program product of claim 6, wherein the 
computer readable program code for low pass filtering further 
comprises computer readable program code for partially 
reconstructing said representation wherein said partially 
reconstructed representation has the form 


Xit it) = Y j c j + r n , 
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where k identifies a cut off frequency, n is the number of 
EMD components, c ■■■ represents IMFs, r„ is the final 
residue and X tt (t) is the extracted trend. 

8. The computer program product of claim 6, further com- 
prising computer readable program code for removing 
extracted said trends from said representation. 

9. The computer program product of claim 8, wherein the 
computer readable program code for removing said trends 
partially reconstructs said representation, said partially 
reconstructed representation has the fomi 


k 

%hk (0 = ^ Cj* 

I 

5 

where k identifies a cut off frequency, c represents IMFs 
and X ;ji (f) is the partially reconstructed representation. 

10. The computer program product of claim 6, further 
to comprising computer readable program code for displaying 
results from said computer readable program code for selec- 
tively filtering IMFs. 



